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ABSTRACT 

We study the long-wave, modulational, stability of steady periodic solutions of the 
Kuramoto-Sivashinsky equation. The analysis is fully nonlinear at first, and can in prin- 
ciple be carried out to all orders in the small parameter, which is the ratio of the spatial 
period to a characteristic length of the envelope perturbations. In the linearized regime we 
recover a high-order version of the results of Frisch, She and Thual, [1], which shows that 
the periodic waves are much more stable than previously expected. 
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1 Introduction 


The Kuramoto Sivashinsky (KS) equation 

tit 4“ XLXLx “b Uxx ”1” WU*xxxx — 0, 

(x,t) G R 1 x R + , (1.1) 

u(x,0) = u 0 {x), 

with v a positive parameter, arises as an amplitude equation in long-wave, weakly nonlinear 
stability analysis, in a variety of applications. It arises, for example, in free surface flow of 
viscous liquids [2, 3, 4, 5], in concentration waves in chemically reacting systems [6, 7, 8], in 
ion diffusion in plasmas [9], in flame propagation [10, 11] and in the dynamics of interfaces 
in two-phase flows [12, 13, 14]. It is the simplest equation with a convective nonlinearity 
and a band of unstable linear modes and therefore a good example on which to apply the 
general notions of inertial manifold theory. This means that the long time dynamics of KS is 
captured well by a finite dimensional dynamical system whose number of degrees of freedom 
is at least as large as the number of linearly unstable modes [15]- [19]. For 2% periodic 
solutions the number of linearly unstable modes is of order i/ -1 / 2 while the best estimate for 
the dimension of the inertial manifold is const • z/~ 21 / 40 , [16], which is quite good except that 
the constant in the estimate is large. 

The KS equation has been studied numerically by many authors [20]-[25]. The main 
objective of the numerical studies is to observe the appearance of low-dimensional dynamic 
behavior, including complicated time-periodic and chaotic behavior. Feigenbaum, [26, 27], 
and others [28, 29], have shown that period doubling and transition to chaos in the iteration 
of one dimensional maps have a universal behavior. This universal behavior is expected 
also in bifurcation of periodic solutions in low-dimensional dynamical systems, for example 
the Lorenz system [30]. Numerical computations [31] have verified this Feigenbaum period 
doubling route to chaos for this system. It is reasonable to expect then that partial differential 
equations whose large time behavior is captured by finite dimensional dynamical systems 
should also exhibit universal period doubling behavior. This was first seen numerically for 
the Ginzburg-Landau equation in [32], for thermal convection in [33] and for the KS equation 
in [34, 35]. Other aspects of the inertial manifold as its dimension increases are explored 
numerically in [36, 37]. 

When v > 1 there are no nontrivial steady, 2n periodic solutions of KS. For v just below 
one there exist steady periodic solutions [38]. Steady solutions with Dirichlet boundary 
conditions and their stability were studied by Novick-Cohen [40]. In this paper we will 
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analyze the modulational stability of steady 2 tt periodic solutions U(y, v) normalized to 
have mean zero 

UUy + Uyy + VlJyyyy = 0, (l.‘2) 

u = r U(y)dy = 0 , U{y) = U(y + 2i r). 

Ztt jo 

For any constants p and c the KS equation (1.1) has the two-parameter family of solutions 

u(t,x ) = pU(p( x — ct)\p 2 v) + c. (1.3) 


which are due to scale and Clallilean invariance of the KS equation. By modulational stability 
we mean the construction and analysis of solutions of KS in all of R near this family, when 
p and c are slowly varying with respect to the period of the steady solution (1.2). In [1], 
the linear modulational stability of U is analyzed by homogenization methods [41]. If t is 
the ratio of the period 2ir to the (long) scale of spatial modulations we must then study the 
behavior of solutions to 


\ rr/ x 2 

rq 4“ U ( )w-p 4~ iixx 4~ e wuxxxx = b 


(1.4) 


as e tends to zero, which is a homogenization problem. This analysis was carried up to second 
order in e in [1] and yields an effective diffusion equation with u dependent diffusion coefficient 
which may be positive, corresponding to stable modulations, or negative, corresponding to 
unstable modulations. However, the stabilizing term u xxxx does not influence the diffusion 
coefficient since it comes in at a higher order in the expansion. It is not surprising therefore 
that the stable range of v obtained in [1] is very small. 

In this paper we analyze the nonlinear modulational stability of U to sufficiently high 
order that the stabilizing term u xxxx is accounted for fully. We use a method similar to the 
one in [42]. In addition to the expansion method, our main result is the determination of a 
much larger range of v than that found in [1] for which there is stability. 


2 Formulation. 

We will construct modulated solutions of KS based on the family of solutions (1.3). Let 

u(t,x) = -U(-‘,p 2 v) + - + w(-,t,x). (2.1) 

e e e e 

where p(t,x), c(t,x), 9 x (t,x) — p(t,x ) and w(y,x,t), which is periodic in the y variable and 
has mean zero, are to be such that for t sufficiently small u solves the suitably scaled KS 
equation 

1 2 
-u t + uu x + u xx -f e vu xxxx = 0, 


( 2 . 2 ) 



In (2.1), e is the ratio of the period ‘ 2 n to a characteristic scale of the modulations. We will 
consider (2.1) at first as a change of variables from u to the triple p, c, w (with 0 X = p and 
the average of w over y, w = 0) and find the exact equations that p , c and w must satisfy 
so that u is a solution of KS. Then we solve the equations for p , c and w approximately for 
small e. The transformed KS equation (2.1) has the form 

p 2 L{w) + F {c) (c,p,w) = 0 (2.3) 


where is a function that depends on p , c and w and their derivatives (it is given in 
Appendix A) and L is the linearized operator for (1.2) with v replaced by p 2 u 

L{(f>) = {U<!>)y + <f>yy + (2.4) 


governing the stability of the cellular solutions U = U(y,p 2 u). As explained further in 
Appendices A and B, in order to be able to solve (2.3) and have a solution w that averages 
to zero in the y variable it is necessary that F ^ have mean zero and that the coefficient of 
U y in F ^ be zero. This leads to two equations for p and c which are 


Pt + ( cp) x + 3e/>xx + e 3 *' 


(3 (pp x )x + PP xx)x + PPx 


= 0 , 


(2.5) 


and 

ct + cc x + U 2 pp x + 2 vUU„p z p x + eU(pw) x + ‘ 2 a/U l/ wp 2 p x 

+ec xx + e 2 ww^ + e 3 vc xxxx = 0, (2.6) 

where the overbar stands for average over the y variable. When these equations are satisfied 
then (2.3) takes the form 


p 2 L(w) + F 0 (c, />; //) + eFi(c, p, u>; ;/) + e 2 F 2 (c, p, w, u) 
+e 3 F 3 (c, p,w, //) + e 4 vw xxxx = 0. 


(2.7) 


where F 0 , F\ and F 3 are given in Appendix A. 

We will also use below the linearized operator L when p — 1 and we will then denote it 

by 

= W)* "b $yy U( i > yyyy (^‘3) 

Some properties of the linearized operator C (or L) are discussed in Appendix B. 
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3 Nonlinear modulation theory. 


The equations (2.5)-(‘2.7) are a coupled nonlinear system equivalent to the KS equation. The 
two modulation equations for p and c are, of course, not closed since they depend on the 
fluctuation term w which in turn depends on p and c. To get a modulation theory we must 
close the p and c equations. We do this by expanding the fluctuation w in a power series in 
e with e — * 0 

w — w 0 + ew i + e 2 w 2 + (3.1) 

The main point of our modulation analysis is that the modulation equations (2.5) and (2.6) 
are exact to all orders in e and we can make the error be formally smaller than the terms 
kept if we retain only the three leading order terms in the expansion of w shown above. 
Substituting (3.1) into (2.5)-(2.7) yields the following equations for the first three terms in 
the expansion of w : 

p 2 L(w 0 ) + F 0 (c, p; v) = 0 (3.2a) 

p 2 L(w i) + Fi(c, p, wo] o) = 0 , (3.26) 


p 2 L(w 2 ) + F 2 (c,p,w 0 ] v) + (( pw\) x U - ( pw x ) x U ) 

+2vp 2 p x (wU u - wU„) + w u + (cwi x ) x + p(w 0 wi y + wiw 0y ) - 2p x wi y 

+ 2pW lxy + 6l/p 2 p X Wlyyy + 4 W 1 x yyy = 0 . (3.2c) 

The functions Fq, F\, F 2 are given in Appendix A. Equations (3.2a-c.) appear to be linear but 
they are in fact nonlinear because they are coupled to the modulation equations (2.5) and 
(2.6). Very little can be said about them at this stage. They must be simplified further. 
In the next two sections we will study the linearized modulation equations and then their 
weakly nonlinear form. 

4 Linear stability theory. 

We now linearize equations (2.5), (2.6) and (3.2a-c) about the exact solution 

p(t,x)=\ , c(f,a;) = 0 , w = 0 . 

while regarding e as a fixed parameter. We obtain the following linear equations for the 
small perturbations to p, c and w 

pt + c x + 3 ep xx + 5e z vp xxxx = 0, (4.1) 
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(4.2) 


c t + (U 2 — 2 uUU u )p x + €w x U + ec xx “h £ VCxxxx — 0} 


C(w 0 ) - 2uU u c x + p x (U 2 - U 2 + 2 v{UU„ - UU V ) 

-\-AvUvy + 8V 2 U l/ yyy + lOuUyyy) = 0, (4.3<l) 

C( Wl) + p xX (-2U - 4 VU„ + 10 vUyy + \2V 2 Uvyy) 

+w ot + ( w 0x U - w 0x U) + 2 w 0xy + Auw 0xyy y = 0 , (4.36) 

C(w 2 ) + w\t + w 0xx + (w\ x U - w Xx U) + 2w x xy 

“f~4 WWlxyyy "f" 6l^ W() XX yy “(" 8lS PxXX^ Vy 0 • (4.3c) 

The solutions of (4.3a-c) can be written in the form 

wo = PxXi(y) + c x X 2 (y ), , (4.4) 

wi = pxxX\\y) + c xx x2 ] {y), (4.5) 

w 2 = pxxxX?\y ) + c xxx x 2 \y ) , (4.6) 

where the functions Xi etc. are functions of y alone and are obtained from 

C(xi) + U 2 — U 2 + 2u{UU l/ — UU V ) + AvUvy + %v 2 U V yy y + 10 vUyyy = 0, (4.7a) 

C(x 2) - 2 vU v = 0, (4.76) 

AxS 1 ’) + (xiU - W) - {TP+ 2vVK)x* + 2 X 1, 

+4^X1 yyy ~ 2C/ — 4 vU u + 10 uU yy + Ylv 2 U vyy = 0 , (4.7c) 

£(Xl^) + (X2 u - X 2 U) — Xl + 2X2 y + ^X2yyy = 0 , (4.7(f) 


£(x? } ) + (xS 1 ^ - xS 1} ) - (■ u 2 + 2 uUU„)xP + 2xff + W$ yy 
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Gi'Xtm ~ 2xi - ("xi txi + 8 v 2 U n = 0, 


(4.7e) 


C(xf ) + (X?U - xPu) - x!" + 2x?„’ 

+4 "xS»v + ( »’X2„ - (t>Xi)Xt = 0 . (4.7/) 

It is convenient to define the weighted averages of the x’s with weight the cellular solution 

U(y). We let 

CVi = Ux 1 , ol -2 = Ux-2 , oc 3 = Ux? , 

«4 = ^X2 ]) 7 «5 = ^Xl 2) , «6 = UX2 ] • (4.8) 

Using (4.4)-(4.6) and the definitions (4.8) in the linearized modulation equations (4.1) and 
(4.2) gives 

pt 4 “ c x + 3 tp xx + 5 t 3 vp X xxx — 0 ? ( 4 * 9 ) 


Cf + (f/ 2 + *2vUUj)p x + COC\p xx + 6(1+ & 2 )c xx 

+ 6 (c^3 /^ xxx + Ot4C xxx ') + 6 Ot§p xxxx + 6 (cTg + V^C XXXX — 0. (4.10) 

Equations (4.9) and (4.10) contain the linear modulational or long- wave stability characteris- 
tics of the cellular steady states U of the Kuramoto-Sivashinsky equation up to order e 3 . The 
parameters of the problem are the coefficients c*i,..., 6V6 which are determined numerically 
as described in Section 6. 

To analyze the stability of U it is sufficient to look for plane-wave solutions of the form 

{p{t, x), c(t, x)) = (p, c)e lkx+wt , (4.11) 

where p, c are constants, k is the wavenumber of a particular disturbance and u is the 
corresponding growth rate. Substituting (4.11) into (4.9) and (4.10) gives the dispersion 
relation that determines w as a function of k 

u> 2 + to [— e(4 + a 2 )k 2 + e 3 (6i / + a 6 )/: 4 ] + k 2 A 2 

+e 2 k 4 ( 3 + Sa 2 — a 3 ) — e 4 k 6 (8u + 3a 6 + 5 a 2 v) + 5e e k 8 v(v + a 6 ) = 0 . (4-12) 

Here 

X 2 =Tji + '2 vTJTT u . 

and in deriving this relation we have anticipated (6.6). We note here that A 2 can be either 
positive or negative, a fact that is central to the modulational stability of the cellular solution. 
It is easy to get the exact solution of the quadratic equation for to. The main feature of our 
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analysis here is the inclusion of higher order terms in e. We will focus on the (linearized) 
structural stability of the celular solution of KS. For k large we have 

u> = ±ik\ + eu>i k 2 -f- . . . , 

where 

u> 1 = i(4 + a 2 ) . (4.13) 

The parameter wj plays the role of a diffusion coefficient. Clearly the system is stable if both 
A 2 > 0 and u\ < 0. Such a two-term result obtained in [1] in their related study. Based on 
the two-term result, they determined a range of values of v for which the cellular solutions 
are linearly stable. It turns out, however, that the extent of this region is very small and 
consequently of limited significance (see [1] and our Section 7). The two-term analysis gives 
an incomplete answer to the linear and weakly nonlinear stability of the system because it 
does not account for the effects of the stabilizing fourth derivative terms in the modulation 
equations. These contributions are easily computed here. We find that there is a correction 
term of 0(ie 2 k 3 ) to u which is dispersive and does not contribute to growth or decay. The 
physically significant for stability correction is obtainable from a large k analysis of the 
dispersion relation (4.13), and is given by 

w = e 3 k 4 ix> 3 , o> 3 = — v — cv 6 or — 5u . (4.14) 

The results (4.13), (4.14) can be used to determine the structural stability of the cellular 
solutions. First, a range of v must be found which satisfies A 2 > 0 so that stable solutions are 
obtained to leading order. The next correction to the real part of u is eh 2 ui \ . For stability 
this term must be negative, and so the range of u for which A 2 > 0 and — u\ > 0 is a 
linear stability region. As mentioned earlier, this range of v is very small because the next 
correction to the real part of u) given by (4.14) is neglected. We have structural stability 
when v is such that — v + «6 > 0. When u>\ is positive but u>z is negative then we have 
structural stability and a finite band of linearly unstable wave numbers. This is a regime of 
v values for which the weakly nonlinear modulation theory of the next section is applicable, 
and it is much larger than the one for which oj\ is negative [1]. The computations leading to 
these conclusions are described in Section 7. 

5 Weakly nonlinear theory. 

To get a suitable weakly nonlinear modulation theory we must return to (2.5), (2.6) and 
(3.2a-c.) and simplify them in a way consistent with the expansion (3.1) without, however, 
altering the construction of the x functions of the linearized theory. This can be achieved by 
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making the size of the perturbations in p, c and w about 1, 0 and 0, respectively, of order 
e 3 , in which case nonlinear terms enter into the modulation equations for c and p and are 
comparable with the stabilizing fourth derivative terms. We therefore let 

p(t,x) = 1 +€ :i p(t,x) , c(t,x) = e 3 c(f,x) , w(t,x) = e 3 t v(t,x) (5.1) 

and substitute this into (2.4)-(2.6) to get the following equations for p and c (with the tildes 
omitted) 

Pt + C x + 3ep xx + e 3 (cp) x + 5 e 3 vp xxxx = 0 , (5.2) 

Ct + A 2 p x + c(l + Of2 )c xx + € 2 GL3 p xxx 

-\-£ 3 (cC x + A 2 ppx) + + v) c xxxx — 0 . (5.3) 

These are the weakly nonlinear modulation equations. The difference between the system 
(5.2) and (5.3) and the linear system (4.10) and (4.11) is the appearance of the nonlinear 
terms of 0(e 3 ) in the former. The constants A 2 , a 2 , a 3 , a 6 for (5.2) and (5.3) are identical 
to the constants of the linear theory in Section 4. 

6 Computation of stability parameters. 

We must solve problems of the form 

£(x) = f{y) , x(y + 2 ?r) = x(y) , x = o , (6.1) 

with f(y) a known 27r-periodic function of mean zero and the operator C defined by (2.8). 
Since x{y) an d /( y) are periodic in y and have mean zero, we can find periodic functions 
?/>(y) and F(y ) such that 

X = Vv , / = Fy . (6.2) 

Substituting into (6.1) and integrating with respect to y gives 

4" tAyj/ 4" V^Pyyyy = F > F — 0, (^*3) 

with ip determined only up to a constant which we take to be zero. It is convenient to 
introduce the operator M. defined by 

■Ad (ip) = Ulpy Ulpy 4" tAyy 4" ^tAyyyy 

Integrating each of the equations (4.7a-f) gives the following sequence of problems 

M{tl> i) 4- J\U 2 - JP)dy 4- 2// J V (UU „ - m)dy 
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+ \OuUyy + 8// 2 U vyy + 4 VU V = 0 , 


(6.3a) 


A4(0- 2 ) — 2v J U„dy = 0 , 


(6.36) 


Mty S !) ) + ftyxyU - tf> ly U)dy - (IP + 2vUU v )xfo + 20 ly + 4 i/0 lwy 


-2 J Udy -Av J U u dy + lOvU y + I2i / 2 U uy = 0 , (6.3c) 

A4(?/>2^) + J (i>2yU - ljj 2 yU)dy — 0! + 20 2j/ + 4 Vljj2yyy = 0 , (6.3d) 

M(j > , (2) ) + /"(^f/ - 0i%)dy - (LP + 2vUU~ L> )^ ) + 2t/> 1 ( * ) 

+ 4// V , iyiv + 6,/ ^iw + 8;/ 2 t/, - 2V>i = 0 , (6.3e) 

+ / y (02y } f/ - 0&^My - 0? } + 20$ 

+4^2^ + 6|/02 W - (t/02y)02 = 0 . (6.3/) 

The function U u is needed in the calculations and is computed by solving 

•Ad(V^) = Uyyyy , U u = fy . (6. 3^) 


We now solve numerically the integrated equations (6.3a-f). We have also solved (5.7a-f) 
by a similar method and with almost identical results but (6.3a-f) is more convenient. The 
central element in the construction of numerical solutions of (6.3a-f) (or (6.1)) is the steady 
state, periodic solutions U of the Kuramoto-Sivashinsky equation, which can be described 
by a very small number of Fourier modes. Since we are interested in a range of values of v 
(approximately 0.3 < v < 1.0) which is above any bifurcations, we know that steady states 
are unique global attractors. The numerical calculations of [34, 35] show that a practical 
estimate of the number of modes required is a few more than i/ -1 / 2 , the number of linearly 
unstable ones. We therefore write 

N 

U(y,v) = Yj U n sm(ny) , 
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where N and U n depend on v and are found numerically. With this approximate solution, 
the linear operator (6.3) becomes a matrix equation for the constants a n ,b n where 

N 

Hy) = LK sin (ny) + b n c.os(ny)). (6.4) 

71=1 

Expanding also the forcing F(y) 

N 

F(y) = E(/» sin (ny) + g n cos (ny)), 

71=1 

the system of linear equations to be solved is 
7u — 1 | N—m y 

E d a iUm-j + E - (i + mfa+mUj) 

i=i 2 i=i z 

+m 2 (i/m 2 - l)a,„ = /„, , (6.5a) 

77i—l i N—m i 

E - E ^ m )bj+mUj) 

j~ 1 j=l 

+m 2 (i/m 2 - l)6 m = y ro , (6.56) 

for l < in < N where the first or second sums are omitted when m = 1 or N respectively. We 
see from (6.5a,b) that the equations for the sine and cosine coefficients of ip are decoupled. In 
particular if /„, = 0 (g m = 0), for m = 1, iV then a m = 0 (b m = 0), m = 1, N also. This holds 
as long as the determinant of the n x n matrix defined by (6.5a) or (6.5b) respectively, is 
non-vanishing. This condition is monitored throughout the computation and is also checked 
analytically near the bifurcation point u = 1 (see Appendix C). The symmetry is due to the 
odd parity of the steady state U(y, v). We expect, therefore, that 

c*i = 0 , c*4 = 0 , <*5 = 0 . (6-6) 

For a given value of v the coefficients U n were obtained numerically by solving the unsteady 
KS equation until a steady state is reached (see [35] for a description of the numerical 
scheme). For values of v near 1, the integration requires longer times to reach a steady state. 
The results reported here have the same accuracy for all values of u. 

With the Fourier coefficients of U ( y ; //) known, the solution vector (aj, a 2 , ..., a n, &i, b n) 
is found by a single matrix inversion for each value of u. The stability constants or,-, i = 1, ..., 6 
are determined easily from their definitions (5.8). The symmetry result (6.6) is confirmed 
by our computations which do not explicitly assume the odd-parity property of U(y) and 
can therefore be extended to general steady states. We note, however, that for the range of 
v studied here, any steady state can be described by the odd-parity solutions by means of a 
suitable translation or Galilean transformation, if necessary. 
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7 Numerical results for linear stability. 

We now present the results of our computations of the stability constants. «i, First 

we consider the behavior of the solutions near the bifurcation point // = 1. This analysis is 
given in Appendix C and the main results are that 

A 2 = -24 + 0(0 , a 2 = -y + 28-16£ (7.1) 

where £ = 1 — v. Analogous results for the asymptotic form of t/> 2 are given in Appendix 
G. Table 1 below compares the asymptotic results for c* 2 with direct computations when 
N =10. 


1 

II 

Exact (N=10) 

Asymptotic 

0.11 

-84.61 

-82.85 

r ‘0.12 

-72.62 

-73.92 

0.13 

-62.84 

-66.39 

0.14 

-54.73 

-59.95 

0.15 

-47.90 

-54.40 


The asymptotic results show that for small £ the solutions are unstable to leading order. 
The next order correction, however, is stable since u>\ < 0 for small 0 For general values of 
v the coefficients are computed numerically as explained in Section 6. In Table 2 below we 
present the results for the parameters A 2 and u\ for a range of values of u. The table also 
gives the values of U. These results were produced with N = 6 or 10 and no change in the 
accuracy was observed. 


V 

A 2 

-«i 

U 

0.89 

-14.59 

40.30 

2.11 

0.85 

-10.88 

21.95 

2.79 

0.80 

-6.84 

10.77 

3.51 

0.75 

-3.33 

4.73 

4.09 

0.70 

-0.287 

1.07 

4.51 

0.65 

2.30 

-1.29 

4.79 

0.60 

4.43 

-2.86 

4.93 

0.55 

6.14 

-3.94 

4.92 

0.50 

7.41 

-4.69 

4.74 

0.45 

8.23 

-5.22 

4.41 

0.40 

8.70 

-5.64 

3.92 

0.35 

8.77 

-6.09 

3.25 

0.30 

8.57 

-6.97 

2.37 
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From Table 2 we see that there is a small region, approximately 0.65 < v < 0.7, where 
both A 2 and —ui\ are positive and so the solutions are stable. The results of Table 2 are 
shown in Figure 1. The stability window (the same as in [1]) is present but is very small. 

When w 3 < 0 then there is a stabilization by the fourth derivative and the structural 
stability window is enlarged to 0.36 < v < 0.7, approximately. Note that — u> 3 becomes 
negative somewhere in 0.36 < v < 0.365 and in 0.695 < v < 0.7, respectively. 


8 Numerical solution of the weakly nonlinear modu- 
lation equations 

We will solve numerically equations (5.2) and (5.3). We consider the following initial value 
problem 

pt + c x + 3ep xx + e 3 (pc) x + ht 3 vp xxxx = 0 , (8.1) 

c t + A 2 p x + e(l + a 2 )c x x + e 2 a 3 p XXX 

+e 3 (cc x + A 2 ppx) + e 3 (o!6 + v)c xxxx = 0 , (8-2) 

p(x,0) = p 0 {x) , c(x,0) = co(x) , {IC) 

where the initial conditions satisfy 

p — > 1 , c — ■> 0 (is |x| — > oo . (7?C) 

We want to see how a localized perturbation in p and c will evolve, when v and e take 
different values. 

For the numerical work it is convenient to transform to new variables p — > 1 + Ap and 
c — » Ac which makes the boundary conditions homogeneous, so that (8.1) and (8.2) become 

Pt ~ b (1 £ 3 ) c x "b Scp X x ~b ^ A[pc) x T 5e vp XXX x — 0 ? (®*^) 

c t + A 2 (l + e 3 )p x + e(l + ct 2 )c x x + e 2 a 3 p xxx 

+c 3 A(cc x + A 2 pp x ) + c 3 (a6 + v)c xxxx = 0 . (8.4) 

The constant A is a measure of the nonlinearity and in the numerical experiments it is such 
that e 3 A = 0.1 or 0.5. 

We will use a spectral method to represent spatial derivatives and a time splitting scheme. 
We will therefore be solving (8.3) and (8.4) on a large (compared to the support of the initial 
data) but finite spatial region, with periodic boundary conditions. We will integrate the 
equations up to the time that p and c disturbances reach the ends of the spatial region, since 
unphysical effects due to the boundary conditions will then be felt which are inconsistent 
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with the modulation theory. For the time splitting scheme, which is second order accurate, 
the nonlinear and linear terms are treated separately. The split equations are 

Pt + e 3 A(pc) x = 0 , (8.5) 

c t + \^ A ( c2 + A'V 2 )x = 0 , (8.6) 

and 

Pt + (1 + ( 3 )c x + 3 ep xx + 5 c 3 ppxxxx = 0 , (8.7) 

Ct + A 2 (l + e 3 )p x + e(l + a 2 )c xx + ( 2 ot 3 p xxx + e 3 (a 6 + v)c xxxx = 0 . (8.8) 

The nonlinear equations (8.5) and (8.6) are a pair of conservation laws. It is easy to see that 
the system is strictly hyperbolic with eigenvalues 

p = e 3 A(c ± A p) , 

the strict hyperbolicity being a consequence of A 2 > 0 for all u of interest. The time integra- 
tion is done by a second order accurate Lax-Wendroff scheme with no artificial viscosity. The 
linear part of the equations is solved exactly in Fourier space. Denoting Fourier transforms 
by a caret, the solutions of (8.7), (8.8) subject to initial conditions (IC) are given by 


p{k,t) = , (8.9) 

c(M ) = , (8.10) 

where At, A 2 , B\, B 2 are constants determined from the initial conditions 

p(k,0) = p 0 {k) , c(k, 0) = co(k) , (8.11) 

Pt(k,0) = (3ek 2 - 5e 3 vk 4 )p 0 (k) - ik(l + e 3 )co(A;) , (8.12) 


c t (k, 0) = (ie 2 k 3 a 3 - ik\ 2 (l + e 3 ))p 0 (k) + (eJfc 3 (l + a 2 ) - e 3 k\u + a 6 ))co(k) , (8.13) 

and given by 

A _Pt—uzpa A Pt — w\po 

At — , A 2 = , 

U)\ — U) 2 U>2 — Wi 

Ct-UzCo Ct-UtCo 

U)\ — lt>2 ^2 — ^1 

The eigenvalues u >\ , u> 2 are the two roots of the dispersion relation 

to 2 + u [-e(4 + a 2 )k 2 + e 3 k 4 (§v + a 6 )] + (1 + e 3 ) 2 k 2 A 2 
+e 2 k 4 (3 + 3a a - (1 + e 3 )cv 3 ) - t 4 k 6 (8v + 3a 6 + 5va 2 ) + 5 t 6 k*v(v + a 6 ) = 0 . (8.14) 
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The solution of the linear part of the problem can then be advanced to the level t + At once 
the spectrum of p and c are known at level t. A further accuracy requirement, particularly 
in the case of imaginary values of u>\ and u> 2 , is that the exponents u>iA t and cu 2 At must 
be smaller in absolute value than a certain tolerance, typically no larger than 0.1 in our 
computations. 

Representative computations are given for the values v = 0.69 and // = 0.55. The 
corresponding values of e used are 0.008 and 0.01 respectively. The linear spectra for these 
parameters computed from the dispersion relation (8.14) are given in Figures 2, 3 and 4 and 
the nonlinear evolutions in Figures 5, 6a-b and 7a-b. 

9 Conclusions 

We have considered the modulational stability of a class of steady spatially periodic solu- 
tions of the Kuramoto-Sivashinsky equation. The theory is initially developed in a fully 
nonlinear framework and analytical and computational results are presented for linear and 
weakly nonlinear regimes. We find that the periodic waves are much more stable than 
previously thought due to a high-order stabilizing correction which follows easily from our 
theory. Nonlinear evolution solutions are also presented and indicate spatial spreading of 
the disturbances. 


Appendices 

A Derivation of the exact nonlinear stability equa- 
tions 


In this Appendix we derive equations (2.5), (2.6) and (3.2a-c). Substituting the exact trans- 
formation (2.1) into the KS equation (2.2) and noting that 

a d 1 . 9 . d 

m~ > m + l% + 2 ' /pp, di;' 

a a i a „ a 

Tx^ Tt + 7% + 2vpp ’Tv' 

with the ^ derivative acting only on (/, yields the transformed KS equation in which e is 
an order parameter. 


~ [(/^i "b p 2 c + 3 Cpp x + t^V^ppxxx + ^r)) U y 
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-\-p ( UUy + Uyy 4" P ^Fyyy /y ) + j P {{^ W )y “I" ^VV P VW HyMy) 
+ (Pt + ( Cp) x )U + (Ot + Cp)Wy + 2Up 2 ( P t + C Px )U v + Avp^PxUvy + %V 2 p 5 p x U vyyy 
+ PPxU 2 + 2 Vp Z p x UU U + 10 Vp 3 p x Uyyy + <H + CC X + f( pw) x U 
+2 eup 2 p x ll u w + ec xx + e 2 ww x + e z vc xxxx ^ 

+ j w t + (cw) x + pwwy + p xx U + 2v{{p 2 p x ) x + ppl)U v + 4uYplU„ 
+p x Wy + 2 pw xy + v(pF + 6{p 2 p x )x)Uyy + 2l V 2 (pF uy + 3 (p 4 px)x + 6p 3 px)Uvyy 

+ 24lS 3 p 5 plUvvyy + §Vp 2 p X Wyyy + 4l//) 3 «; *3/3/3/] 

to., + 2 v\pF v + F? + pp x F)U vy + \v z ( P F vv + 3 (p 4 p 2 x ) x + PPxF^U^ 

+ 32 U 4 p 5 p 3 x U l/l/ ^y + \2vpp x w x yy + 6 P f> 2 W xx yy + // F Wyy 

+ev [p xxxx U + 2 v(F x + pp x p xxx )u u + 4 v 2 {F^ + pp x F u )U vv 
+8 v 3 ((p 4 pl) x + ppxF vv )U uvv + 16 u 4 p 5 p 4 x U^ “|“ ^pWxxxy 

~V^Pxx^xy + ^ Px^xxy “I” Pxxx^y\ ^ WU) xxxx 0 * 


(41) 

Here 


F 3 [ppx)x + PPxX 1 

F — (/? Px)xx + ( PPx)x + PPxpxx i 
F uy — 2 (p 3 p x ) x + 4 p 2 p 2 + p(p 2 p x )x > 

F 1 ' 1 ' = ( p 3 pi)x + ppx{p 2 px)x + P*pl • 


Equation (Al) above lias been rearranged in ascending powers of e and certain terms 
have been grouped together as we explain next. The objective is to derive exact nonlinear 
modulational equations which hold in the asymptotic regime e <C 1 without the need of 
cumbersome solvability conditions at each order of the expansion. The modulation equations 
act as all order solvability conditions. This is achieved by making the fluctuation w i have 
zero mean. Since the determination of w c depends on solutions of an equation of the form 
L(w c ) = f(y ) where f(y) is a known forcing function, it is necessary that / have mean zero 
and that w ( has zero mean if f(y) does not contain terms proportional to lJ y . All terms 
proportional to U y , therefore, are grouped together, as shown in (Al) above, and c , p are 
chosen so that the mean of / and the coefficient of U y vanishes for all x and t. The equations 


, , x . .3 (3 {pp x )x + pp xx)x + PPxxx 

pt + \Cp) x + oepxx + € V 
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c t + cc x + U 2 pp x + 2vUU v p 3 p x + eU(pw) x + 2euU l/ wp 2 p x 


~\~€C xx ”1” ^ WU) X “j" € vc xxxx — 0 • (A?) 

Note that (2.3a) is obtained by dividing (A2) by p and differenting with respect to x. With 
the above choice of modulation functions, the equation that governs the fluctuation w be- 
comes 

p 2 L{w) + F 0 + cF\ + c 2 F 2 4* e z vFz + = 0. 

where the functions F 0 , F\ , F ' 2 , Fz are 

Fq — 2 Vp C x LI 1/ 4" 10 Up p x U yyy -f” 8 V p P x F i/yyy 

+Avp 3 p x U l/y + PPx {U' 2 -TP) + 2 up 3 p x {UU u - ulh) , 

Fi = ( pw) x U — ( pw) x U -|- 2up 2 p x (wU l/ - wU„) + w t + ( cw) x 
4 - pwwy — 2 p x w y - dv^pxxU,, — 3 p xx + p xx U 4- 2u(pp 2 x + (p 2 p x )x)U l/ 

-\-Av 2 p 3 p 2 x U w + u(pF + 6(p 2 px)x)Uyy + 2u 2 (pF l>y + 3 (p 4 p x )x + ^p 3 pi)^vyy 
4" 24// P p X U isi/yy 4" 2 ptV X y 4 6//p' Px^yyy 4* 4 Up W X yyy , 


F 2 = 2u 2 (pF 1 ' + F^ + pp.F^y + ^ipF^ + 3(p 4 pl) x + pp.F^U^y 
4-32/z 4 p $ p z x U vvvy 4- vo xx 4- ww x — wwZ 4- \2upp x w xyy 4- §vp 2 w xxyy 4- v{ppxx 4- 3 (pp x ) x )w y y , 


II 

G? 

( PPxxx F x \ 

pxxx x 1 1 

U + 2// 

T?v 1 2 ( PPxxx ~b Fx\ 

F x “t" ppxpxxx P 1 1 

1 

V p Jx\ 


L \ P Jx\ 


IK 


4-4 v 2 (F^ + pp x F v )\J vv 4- 8 v 3 {{p 4 p 3 x ) x + ppxF^U^ 4- 1 6// 4 p 5 ^f7, 


+ 


Pxxx 


PPxxx + F x 


UJy “h ^p^xxxy 4 Pxx^xy "1“ ^ Px^xxy j 


B Properties of the linear operator £. 

At several points in the analysis we have to solve 

C{w) = h, 


(B 1) 
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where, omitting x and t dependence, h = h(y ) is a given forcing function with mean zero 
and w =~w(y) is also a periodic function with mean zero. The operator £ is defined by (2.8). 
It is clear that 

40 = Vy 

and so to make the solution of (Bl) have mean zero we must set to zero the coefficient of U y 
in h. Thus, the solvability condition for (Bl) is 

• h = 0, and 

• The coefficient of U y in h (if any) is zero. 

C Stability analysis near v = 1. 

We are interested in the modulational stability of the cellular solution U for values of v just 
below 1. Let // = 1 — £ where ( <1. A two-term Fourier expansion is sufficient to describe 
the steady state U(y) : 

U (y, 0 = U x sin (y) + U 2 sin(2 y) + . . . , 

Ux = -2 12 1/2 £ 1/2 , U 2 = — 2£ . (CT) 

We also need U u which from (Cl) is 

Uv{y\£) = 12 1/2 ^ _1/2 sin(t/) + 2sm(2j/) + ... . (C2) 

To fix matters, consider the asymptotic solutions for rp 2 from equation (6.3b). This solution 
leads directly to the value of a 2 which is one of the important stability parameters. The 
solution -02 is 

^2 {y) = ax sin(y) + a 2 sin(2y) + &i cos(y) -f b 2 cos(2 y) + ... . (C 3) 

Since the forcing terms of (6.3b) can be expanded in a series of cosines, we expect ax = a 2 = 0. 
This holds for general values of v (verified numerically), and it is shown here that the result 
also holds near v = 1. Substituting (C3) into (6.5a) and (6.5b) with N = 2 and m = 1,2 
yields the decoupled linear equations 


( 2^2 + (4 — l))fli — Uxa 2 — 0 , 

(C4.1) 

If/iai +4(4)/ - l)o 3 = 0 , 

(C4.2) 


17 



(C75.1) 


{-\ih + (u - l))6i - U\ 6 2 = -2(1 - Ol2 1/2 r 1/2 , 

+ 4(4// - 1)6 2 = -2(1 - 0 . ((75.2) 

The determinant of the coefficient matrix of the first system is D = 32£ 2 , and so the trivial 
solutions <ii = 02 = 0 follow as expected. The second system gives 

»■ = - §1 + 18 < ,/! ) > ^ = 4 + 1 ■ (C6) 

These expressions are exact in £. From the definition (4.9) of c *2 and using X 2 = i’-iy we find 

a 2 = -\lhh ~U 2 b2 = -^ + 28- 16{ . (Cl) 

We note that this result can also be obtained from an asymptotic analysis of the differenti- 
ated equation (5.7b) for X 2 - Comparison of the asymptotic expression (C6) with the direct 
simulations shows that they are in excellent agreement (see Section 7). 
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Stability results 



V 


• Figure 1 Linear stability characteristics, e = 0.025. Structural stability is achieved 
when A 2 and — 0 J 3 (from the fourth derivative terms) are positive. If we stop at 0(k 2 ), 
a small structural stability window is obtained when A 2 and — U\ are positive, as seen 
in the figure. Inclusion of higher order terms extends the structural stability window 
significantly. 
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Spectrum of most unstable mode 

vw0.69 c=C.9S 



• Figure 2 The most unstable eigenvalue computed from the exact dispersion relation 
(8.14) or equivalently (4,12) (the difference is a scaling factor due to a change of 
variables). The parameters are u = 0.69, e = 0.95 corresponding to a stability window 
where both the 0{k 2 ) and 0{k 4 ) terms are dissipative. 


Spectrum of most unstable mode 


v»0.69 e=*0.008 



• Figure 3 As in Figure 2, but e = 0.008. For this e a large but finite band of unstable 
modes is obtained. 
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Spectrum of most unstable mode 


v=0_55 €=0.01 



k 

• Figure 4 As in Figure 2, but v = 0.55, e = 0.01. For this v the 0(k 2 ) term is unstable 
and the 0(k 4 ) term is dissipative. There is a small band of stable modes near A; = 0 
and a large band of unstable modes. 


Weakly nonlinear evolution 

v*Q.69 c*0JJ 



• Figure 5 Solution of the nonlinear system (8.3), (8.4). The initial condition for p 
and c is a Gaussian of amplitude 0.25 centered at the midpoint of the domain. The 
evolution of p is shown for v = 0.69, e = 0.95. The corresponding c decays fast to a 
steady state (equal to the mean of the initial condition) by the end of the computation. 
The spectrum corresponding to this run is shown in Figure 2. 
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Figure 7a As in Figure 5, evolution of p for v = 0.55, e = 0.01. The spectrum is 
shown in Figure 4. 


Weakly nonlinear evolution 


v*OJ5 e*0.01 
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